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£SJ , We present Nt — 4 lattice results for the equation of state of 2 + 1 flavour staggered, 

dynamical QCD at finite temperature and chemical potential. We use the overlap improving 
multi-parameter reweighting technique to extend the equation of state for non-vanishing 
chemical potentials. The results are obtained along the line of constant physics. Our physical 
parameters extend in temperature and baryon chemical potential upto ~ 500 — 600 MeV. 



§1. Introduction 



> . 

£N| . QCD at finite T and/or \i is of special importance since it can be used to describe 

the early universe, neutron stars and also heavy ion collisions. Present and future 
heavy ion collisions are carried out at CERN, in Brookhaven and at GSI to detect and 
study experimentally the QGP phase, i.e. QCD at large temperature and moderate 
chemical potentials. 1 )' 3 - 1 It is very important to understand the theoretical grounds 
of the underlying physics. First principle answers can be gained from lattice QCD 
calculations e.g. to obtain the equation of state (EoS). 

The experiments are carried out at /i ^ but unfortunately until recently the 
lattice results were limited to /i = 0. Though lattice QCD can be easily formulated 
at non- vanishing chemical potentials 4 )' 5 ) we cannot use Monte-Carlo simulations at 
/i 7^ as the determinant of the Euclidean Dirac operator and so the functional 
measure becomes complex. 

Recently two of us proposed a new technique, 6 ) the so-called overlap improving 
multi-parameter reweighting method to study lattice QCD at finite [i . This pro- 
cedure proved to be good enough to give the phase boundary on the T-/x plane for 
four flavours, 6 ) for 2+1 flavours 7 ) and the equation of state. 8 )' 9 ) Essentially the same 
technique was used succesfully by other studies ; 10 )~ 12 ) however, instead of evaluating 
the fermionic determinant exactly it can be approximated by its Taylor series with 
respect to /U. Other approaches, like simulations at imaginary chemical potential 
and analytic continuation led to results that are in good agreement with those of our 
method. 13 H 5 ) 

In this paper we determine the EoS on the line of constant physics (LCP). An 
LCP can be defined by a fixed ratio of the strange quark mass (m s ) and the light 
quark masses (m u d) to the [i = transition temperature (T c ). Our parameter choice 
approximately corresponds to the physical strange quark mass. However, the ratio 
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of the pion mass (m n ) and the rho mass (m p ) is around 0.5-0.75, which is roughly 3 
times larger than its physical value. In our lattice analysis we use 2 + 1 flavour QCD 
with dynamical staggered quarks. The determination of the EoS at finite chemical 
potential needs several observables at non-vanishing ^-s. These are produced by the 
use of the multi-parameter reweighting method. We employ the integral method to 
calculate the pressure. 16 ) 

The paper is organized as follows. In Section 2 we summarize the lattice param- 
eters and the technique by which the lines of constant physics can be determined. 
Section 3 presents the equation of state at vanishing chemical potential. Sections 4 
deals with the question how to reweight into the region of /j, + 1 and how to estimate 
the error of the reweighted quantities. In Section 5 we give the equation of state for 
non-vanishing chemical potential and temperature. Those who are not interested in 
the details of the lattice techniques should simply omit Sections 2-4 and jump to 
Section 5, or refer to. 8 ) Finally, Section 6 contains a summary and the conclusions. 
The details of this work can be found in. 9 ) 

§2. Lattice parameters and the line of constant physics 

In this paper we use 2 + 1 flavour dynamical QCD with unimproved staggered 
action. Simulations are done for the equation of state along two different lines of 
constant physics and at 14 different temperatures. The temperature range spans up 
to 3T C . In physical units our parameters correspond to pion to rho mass ratio of 
m n /m p ~ 0.5 — 0.75 and lattice spacings of a ~ 0.12 — 0.35 fm. 

The finite temperature contributions to the EoS are obtained on 4 • 8 3 , 4 • 10 3 
and 4 • 12 3 lattices, which can be used to extrapolate into the thermodynamical 
limit (we usually call them hot lattices). On these lattices we determine not only 
the usual observables (plaquette, Polyakov line, chiral condensates) but also the 
determinant of the fermion matrix and the baryon density (%) at finite [i. 10000 — 
20000 trajectories are simulated at each bare parameter set. Plaquettes, Polyakov 
lines and the chiral condensates are measured at each trajectories whereas the CPU 
demanding determinants and related quantities are evaulated at every 30 trajectories. 
For our parameters the CPU time used for the production of configurations is of the 
same order of magnitude as the CPU time used for calculating the determinants. 

Since we usually move along the line of constant physics by changing the lattice 
spacing a and keeping the masses fixed we will explicitely write out the lattice spacing 
a in our formulas. In this paper we study lattices with isotropic couplings. We write 
HB for the baryonic chemical potential, whereas for the quark chemical potential 
(u, d quarks) we use the notation /x. Similarly, the baryon density is denoted by ns 
and the light quark density by n. 

In the remaining part of this section we discuss the role of LCP when determining 
the EoS in pure gauge theory and in dynamical QCD. After that we determine the 
lines of constant physics, along which our simulations are done. 

In order to detetermine the temperature T = l/[N t a((3)] of the pure gauge 
theory, we have to compute the lattice spacing (a) as a function of the gauge coupling 
(/?). In the d dimensional space of the bare parameters one defines d appropriately 
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chosen quantities. The LCP is given by d — 1 constraints and it is parametrized by a 
non-constrained combination of the above quantities. For the 2 + 1 flavour staggered 
action we have three bare parameters (/3, and two masses, m ud ,m s ). Thus, we need 
two constraints. There are several possibilities for these constraints and consequently 
there are many ways to define an LCP. A convenient choice for two of the three 
quantities can be the bare quark masses (m u d and m s ). A more physical possibility 
is to use the pion and kaon masses (m,,-, mif). 

In our analysis first we use the bare quark masses (777,^ and 777 s ) and the tran- 
sition temperature to define an LCP. In this paper we use two*) LCP's (LCPi and 
LCP2). The conditions 



are taken as the constraints for LCPi and LCP2, respectively. For both LCP's 
we determined four different transition couplings (f3 c ) by susceptibility peaks on 
Nt = l/(T c o) = 4, 6, 8 and 10 lattices with spatial extensions N s > 2>Nt/2 and quark 
masses given by eq. (|2-1|) . The quark masses or the transition gauge couplings can 
be used to parametrize the LCP's. 

By the finite temperature technique, described above, only a few points of the 
LCP's can be obtained. To interpolate (5{a) between these points (and extrapolate 
slightly away from them) we use the renormalization group inspired ansatz proposed 
by Allton. 17 ) A particularly illustrative parametrization is obtained by inverting 
eq. (|2Tj) and using Nt as a continuous parameter. Figure ^ shows LCPi and LCP2 
with our simulation points. The simulation points in the "non-LCP" approach - often 
used in the literature - are also shown. Note that even though the determination 
of the LCPi and LCP2 are done on finite temperature lattices, the obtained bare 
parameters are used in the rest of the paper for T = and simulations. 

§3. Equation of state along the lines of constant physics (LCP) at = 

In previous studies of the EoS with staggered quark actions, the pressure and 
the energy density were determined as functions of the temperature for fixed value 
of the bare quark mass 777 g a in the lattice action. 18 )~ 20 ) In these studies a fixed 
N was used (e.g. Nt = 4 or 6) at different temperatures. Since T = l/{N t a) the 
temperature is set by the lattice spacing which changes with (3. This convenient, 
fixed bare m q a choice leads to a system which has larger and larger physical quark 
masses at decreasing lattice spacings (thus, at increasing temperatures). Increasing 
physical quark masses with increasing temperatures could result in systematic errors 
of the EoS. 

Clearly, instead of this sort of analysis (in the rest of the paper we refer to it as 
"non-LCP approach") one intends to study the temperature dependence of a system 
with fixed physical observables, therefore on an LCP. 

*^ As we will see later, two LCP's are needed for the determination of the EoS at finite chemical 
potential. 



m ud = 0.48T C = 0.48/ (N t a) and 
m ud = 0.384T C = 0.384/(iV» and 



777,5 = 2.08 • m ud 
m s = 2.08 • m ud 



(2-1) 
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Fig. 1. The lines of constant physics (LCPi and LCP2) on the j3 vs. m^a plane. The strange 
quark mass is given by m s — 2.08m u[ j for both LCP's. The simulation points are shown by 
squares/triangles and connected by dashed/dotted lines for LCP1/LCP2, respectively. The 
diamonds along a horizontal line represent the simulation points in the "non-LCP" approach. 
Additional 4 diamonds in the vertical direction show the simulation points used to test the path 
independence of the integral method. 



In our analysis we use full QCD with staggered quarks along the LCP and 
compare these results with those of the "non-LCP approach" . 

Now, we briefly review the basic formulas and emphasize the issues related to 
the EoS determination along an LCP. 

The energy density and pressure are defined in terms of the free-energy density 

(/): 



e(T) = f-T^L, p(T) = -f. 



(3-1) 



Expressing the free energy in terms of the partition function (/ = —T/V log Z 
-Td(\ogZ)/dV) we have: 

T 2 dlog Z 



e(T) 



p(T) = T 



d{\o g z) 



V dT iy ' dV 

The temperature and volume are connected to this lattice spacing by 

1 



(3-2) 



aN t 



(3-3) 



Inspecting eas. (|3-21 13-3|1 we see that (e — 3p)/T 4 is directly proportional to the 
total derivative of log Z with respect to the lattice spacing: 

e-3p Nf d(logZ) 



7%4 



da 



(3-4) 
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Here, the derivative with respect to a is defined along the LCP, which means that 
only the lattice spacing changes and the physics (in our case m g /T c ) remains the 
same. We can write: 

d _ d(3 d d{m q a) d 

da da d/3 ^ da d(m q a) ' 

Since the LCP is defined by m q /T c = const., the partial derivative d{m q a) / da be- 
comes simply m q . The derivatives of log Z with respect to (3 and m q are the plaquette 
and tytyq averages multiplied by the lattice volume. We get: 



LZp. = -Nf a (Pl^ 
T 4 * \ 8a 



LCP 



q ) 



(3-6) 



The pressure is usually determined by the integral method. 16 ^ The pressure 
is simply proportional to log Z, however it cannot be measured directly. One can 
determine its partial derivatives with respect to the bare parameters. Thus, we can 
write: 



n[ r^ a) ( diogz/dp 

NiJ { p , mqOa) d[(j > mqa) \ dlogZ/d(m q a) 



(3-7) 



Since the integrand is the gradient of log Z, the result is by definition independent 
of the integration path (we explicitely checked this path independence). For the 
substracted vacuum term we used the zero temperature pressure, i.e. the same 
integral on Nto = N s lattices. The lower limits of the integrations (indicated by 
(3q and m q o) were set sufficiently below the transition point. By this choice the 
pressure becomes independent of the starting point (in other words it vanishes at 
vanishing temperature). In the case of 2+1 staggered QCD eq. (|3-7|) can be rewritten 
appropriately and the pressure is given by 

v r(l3,m q a) ( (PI) 

- -Nf / d(p,m ud a,m s a) (m ud ) I , (3-8) 



T~>4 

J(/3 ,m q0 a) \ (}J/\J/ s 

where we use the following notation for subtracting the vacuum term: 

(0((3, m)) = 0(P, m) T ^ - 0(f3, m) T=0 . (3-9) 

Figure |^1 shows the EoS at vanishing chemical potential on Nt = 4 lattices for 
LCP2 and for the non-LCP approach. The pressure and e — 3p are presented as a 
function of the temperature. The parameters of LCP2 and those of the non-LCP 
approaches coincide at T = T c . The Stefan-Boltzmann limit valid for Nt = 4 lattices 
is also shown. 

It can be seen that the EoS along the LCP and the "non-LCP approach" differ 
from each other at high temperatures. This is an obvious consequence of the fact 
that in the "non-LCP approach" the bare quark mass increases linearly with the 
temperature. The "LCP" pressure is much closer to the SB limit. 
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Fig. 2. The equation of state at fi—0. (a) The left panel shows the pressure p, as a function of 
the temperature. All quantities are normalised by T 4 . In order to lead the eye the solid line 
connects the data points obtained along the LCP2, whereas the dashed line connects the data 
points obtained in the "non-LCP approach" (see text). The Stefan-Boltzmann limit is also 
shown by an arrow for Nt—i lattices. The EoS along the LCP and the "non-LCP approach" 
differ from each other at high T. (b) The same for e-3p. 



§4. Reliability of reweighting and the best reweighting lines 

The aim of this section is to study the reliability of the multi-parameter reweight- 
ing (for another study of reweighting see 21 )). We also determine its region of validity 
by a suitably estimated error. To start with let us briefly review the multi-parameter 
reweighting. 

As proposed in 6 ) one can identically rewrite the partition function in the form: 
Z{m,n,0) = JvUexp[-S bos (p ,U)]detM(m ^ = 0,U) (4-1) 

where U denotes the gauge field links and M is the fermion matrix *' . The chemical 
potential \i is included as exp(o/x) and exp(— a\i) multiplicative factors of the forward 
and backward timelike links, respectively. In this approach we treat the terms in the 
curly bracket as an observable - which is measured on each independent configura- 
tion, and can be interpreted as a weight - and the rest as the measure. Thus the 
simulation can be performed at fx = and at some f% and mo values (Monte-Carlo 
parameter set). By using the reweighting formula (|4Tjl one obtains the partition 
function at another set of parameters, thus at fi 7^ 0, (3 7^ f% or even at m / mo 

*' For nf 7^ 4 staggered dynamical QCD one simply takes fractional powers of the fermion 
determinant. 



The QCD equation of state at finite T/fjL on the lattice 



7 



(target parameter set). *) 

Expectation values of observables can be determined by the above technique. 
In terms of the weights (i.e. the expression in the curly bracket of eq. (|4-1|) ) the 
averages can be determined as: 

}2{w{f3,n,m,U)} 

Now we present an error estimate of the reweighting procedure which also shows 
how far we can reweight in the target parameter space **\ 

The steps of the new procedure are as follows. First, we assign to each config- 
uration of our initial sample the weight w valid at the chosen target parameter set 
i.e. 

exp {AS • V • (PI) + ^pndet M(ji) - lndetM(// = 0)]} (4-3) 



w 



where V is the lattice volume. Next, we carry out Metropolis-like accept /reject steps 
with the series of these new weights. This procedure generates a new, somewhat 
smaller, sample. The configurations of this new sample are taken with a unit weight 
to calculate expectation values, variances and integrated autocorrelation times 27 ) 
(the latter grows at every rejection due to the repetition of certain configurations). 
The expectation values are taken from eq. (|4-2|) . whereas the errors are estimated 
from the new procedure. 

We still have to clarify two important points: thermalization and the possible 
lack of relevant configurations. We solved the first problem by defining a thermal- 
ization segment at the beginning of every newly generated sample which we cut off 
from the sample before calculating the expectation values and the errors. An obvi- 
ous assumption is to claim that the already thermalized sample containing valuable 
information starts with the first different configuration right after the one with the 
largest weight. This ensures that if there is only one configuration in the initial 
sample which "counts" at the target parameters then this information will not be 
lost. The second problem cannot be solved perfectly. This problem occurs e.g. when 
a phase transition is very strong. 

To illustrate the new technique let us take a look at the nf = 4 flavour case at 
m q a = 0.05 bare quark mass on 4 • 6 3 size lattice at imaginary chemical potential. 
Note that for purely imaginary chemical potentials direct simulation is possible, 
therefore it is possible to check the validity of any possible error estimation method. 
We carried out simulations***) at Im(/i) = in the phase transition point, i.e. at 
j3 = 5.04 (« 2500 independent configurations) and at (3 = 5.085 (« 7000 independent 
configurations). From these two starting points with the use of the reweighting 
we tried to predict the plaquette (PI) and the Tpip expectation values and their 
uncertainties at f3 = 5.085 and Im(/x) ^ 0, that is at the target, imaginary chemical 
potential values. We calculated the plaquette and the rpif) expectation values by 



*' Note the appearent similarity of the present method with those of. 22 - ) ~ 26 ' ) 
**' For other techniques to estimate the errors of the reweighting method see e.g. 27 '' 28 ' 
***) This parameter set is identical to the one used in 6 ' 
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Fig. 3. (a) The left panel shows a well functioning error estimate. The squares correspond to 
direct simulations (out of 60000 configurations), the crosses denote the results of the multi- 
parameter reweighting method (out of m 2500 independent configurations) while the circles 
are the points of the Glasgow- type reweighting (out of ps 7000 independent configurations), 
(b) In the right panel the meaning of the symbols are unchanged but the sample sizes for the 
reweighting techniques are decreased to 1200 independent configurations. Noticeably, in case 
of the Glasgow-type reweighting the fact that the sample does not contain any configuration 
from the target phase causes systematic errors in the expectation values and mainly in their 
uncertainties. This is because in case of strong phase transitions - like in our case - it is far 
too difficult to define a reliable thermalization stage by virtue of the weights. By increasing the 
sample size helps in this situation leading to the results of the left panel. 

(|4-2|) at the target points, and we also used the new Metropolis-type method defined 
above which leads to very similar results. 

The results are shown in Figure |3] (see explanation there) . The infinitely large 
errors of the single-parameter reweighting (Glasgow-method) indicate that the whole 
sample is thermalization, that is it does not provide information about the expecta- 
tion value in the required point. In the right panel of Figure |3] the second problem 
(lack of relevant configuration) mentioned above is seen. 

When we determine the EoS the jackknife errors are suitable to estimate the 
uncertainties of the reweighted quantities in an appropriate region. We used the 
new error estimates only to provide the limit of the applicability of the reweighting 
procedure and of the jackknife method. 

We can define reweighting lines on the /3-(i plane so as to make the least possible 
mistake during the reweighting procedure. To do this we introduce the notion of 
overlap measure which we denote by a. The overlap measure is the normalised 
number of different configurations in the sample created with the Metropolis-type 
reweighting after cutting off the thermalization. We plotted the contour lines of a 
in the left panel of Figure |1] The dotted areas are unattainable, that means here 
the overlaps vanish, the errors are infinitely large. The best reweighting line can be 
defined for each simulation point. For a given value of \x we choose (3 so that a be 
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Fig. 4. (a) The left panel shows the real chemical potential-/? plane. 33000 configurations were 
simulated at the parameter set at the critical f3 in the rij — 2 + 1 flavour case. The dotted 
lines are the contours of the constant overlap. The dotted area is the unknown territory where 
the overlap vanishes. The solid line is the phase transition line determined by the peaks of 
susceptibility, (b) In the right panel the volume and the fi dependence of the overlap (a) is 
shown. Upper curves correspond to smaller lattice sizes, 4-6 3 , 4-8 3 , 4-10 3 and 4-12 3 respectively. 



maximal. The points of the best reweighting lines are given by the rightmost points 
of the contours of the overlap in Fig. 0] (a) . 

Then the best reweighting lines are the contours of constant overlap or equiva- 
lently of constant error. 

The right panel of Figure |1] shows the fi dependence of the overlap at fixed (3 
and quark mass parameters for different volumes (V = 4 • 6 3 , 4 • 8 3 , 4 • 10 3 and 
4 • 12 3 ). As expected, for fixed \x larger volumes result in worse overlap. One can 
define the "half-width" (/ii/2) of the /U dependence by the chemical potential value at 
which a = 1/2. One observes an approximate scaling behaviour for the half-width: 
^1/2 with 7 1/3. 

It is obvious that the two-parameter reweighting used previously does not follow 
the LCP (P gets smaller but the quark mass remains rriQa). Nevertheless, the best 
reweighting line along the LCP can be determined by two techniques. One of them 
is the three- parameter reweighting, the other one is the interpolating method. 

As it can be seen in the left panel of Figure HJ the change in (3 is not very large 
for the two-parameter reweighting. Therefore, one can remain on the LCP by a 
simultaneous, small change of the mass parameter of the lattice action. This results 
in a three-parameter reweighting (reweighting in ma, (3 and fia). Similarly to the 
two-parameter reweighting one can construct the best three- parameter reweighting 
line. 

Another possibility to stay on the line of constant physics at finite /U is the 
interpolating technique. One uses the two-parameter reweighting for two LCP's and 
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interpolates between them. The result of this method and the predicition of the 
three-parameter reweighting agree quite well. This indicates that the requirement 
for the best overlap selects the same weight lines even for rather different methods. 

§5. Equation of state at non-vanishing chemical potential 

In this section we study the EoS at finite chemical potential. Since we are 
interested in the physics of finite baryon density we use \x = \i u = fid 7^ for the two 
light quarks and fj, s = for the strange quark. 

The pressure (p) can be obtained from the partition function as p=T-d log Z/dV 
which can be written as p=(T/V) ■ \ogZ for large homogeneous systems. On the 
lattice we can only determine the derivatives of log Z with respect to the parameters 
of the action (J3, m, //), so p can be written as a contour integral: 16 ^ 

The integral is by definition independent of the integration path. The chosen inte- 
gration paths are shown on Fig 03 

The energy density can be written as e = (T 2 /V) ■ d(log Z) / dT + (fiT/V) ■ 
d(log Z)/dfi. By changing the lattice spacing T and V are simultaneously varied. 
The special combination e — 3p contains only derivatives with respect to a and 

e — 3p a dlog(Z) 

T 4 ~ ~ TH? Qa 

We present lattice results on Ap(fi,T) = p(fx 7^ 0, T) — p(fi = 0, T), e(fi,T)- 
3p(n,T) and ns(//, T). Our statistical errorbars are also shown. They are rather 
small, in many cases they are even smaller than the thickness of the lines. 

On the left panel of Fig. Elwe present Ap/T 4 for five different fi values. On 
the right panel normalisation is done by Ap SB , which is Ap(/j,,T — > 00). Notice 
the interesting scaling behaviour. Ap/Ap SB depends only on T and it is practically 
independent from fi in the analysed region. The left panel of Fig. [7| shows e-3p 
normalised by T 4 , which tends to zero for large T. The right panel of Fig. [7| gives 
the dimensionless baryonic density as a function of T/T c for different fj,-s. 

The error coming from reweighting has been discussed previously. Another 
source of error is the finiteness of the physical volume. The volume dependence 
of physical observables is smaller than the statistical errors for the plaquette average 
or quark number density. 

§6. Conclusions, outlook 

We studied the thermodynamical properties of QCD at finite chemical potential. 
We used the overlap improving multi-parameter reweighting method. Our primary 
goal was to determine the equation of state (EoS) on the line of constant physics 
(LCP) at finite temperature and chemical potential. 



+ 



/f gflogg) 
T 3 V dfx 



(5-2) 
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Fig. 5. Illustration of the integral method at finite chemical potential (left panel). The solid lines 
are fi — const lines on the fxa — f3 plane. Dashed lines are the best reweighting lines starting 
from different simulation points. Arrows show the path of integration we used when evaluating 
eq. (7.6). 



We have pointed out that even at /i=0 the EoS depends on the fact whether 
we are on an LCP or not. According to our findings pressure and e - 3p (interacion 
measure) on the LCP have different high temperature behaviour than in the "non- 
LCP approach". 

We discussed the reliability of the reweighting technique. We introduced an 
error estimate, which successfully shows the limits of the method yielding infinite 
errors in the parameter regions, where reweighting gives wrong results. We showed 
how to define and determine the best weight lines on the /J.—/3 plane. 

We discussed the two-parameter reweighting technique. Two techniques were 
presented (three-parameter reweighting and the interpolating method) to stay on 
the LCP even when reweighting to non- vanishing chemical potentials. 

We calculated the thermodynamic equations for fi ^ and determined the EoS 
along an LCP. We presented lattice data on the pressure, the interaction-measure 
and the baryon number density as a function of temperature and chemical potential. 
The physical range of our analysis extended upto 500 — 600 MeV in temperature and 
baryon chemical potential as well. 

Clearly much more work is needed to get the final form of non-perturbative EoS 
of QCD. Extrapolation to the thermodynamic and continuum limits is a very CPU 
demanding task in the /i ^ case. Physical m n /m p ratio should be reached by 
decreasing the light quark mass. Finally, renormalised LCP's should be used when 
evaluating thermodynamic quantites. 
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Fig. 6. The equation of state at, = 100, 210, 330, 410 and 530 MeV. The left panel shows 
the pressure difference between the jj, = and fj, ^ cases normalized by T 4 , whereas on the 
right panel the normalization is done by the iVt = 4 lattice Stefan-Boltzman limit. Note, that 
Ap/Ap SB seems to show some scaling behaviour (it depends mostly on the temperature; whereas 
its dependence on \x is much weaker). Thus, the fi dependence of Ap is almost completely given 
by the n dependence of the free gas. 
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Fig. 7. (a) (e — 3p)/T 4 (b) Dimensionless baryon number density as a function of T/T c at 
Ats=100,210,330,410 MeV and ^ s =530 MeV. 
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